##################################################################################################################################################################################################################
# Replication file 2 for Democrats versus Democracy: The Southern Response to the Voting Rights Act
# This file produces all of the main analysis for the paper and SI. Run this after you run gen_county_data.R
##################################################################################################################################################################################################################

# Set Filepath to the Directory That Contains this Code
  setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

##################################################################################################################################################################################################################
# Load Libraries, Clear Workspace
##################################################################################################################################################################################################################
# Formatting
  library(labelled)
  library(dplyr)
  library(tidyverse)
  library(haven)
  library(stringr)

# Plotting & Tables
  library(ggplot2)
  library(stargazer)
  library(texreg)
  library(RColorBrewer)
  library(data.table)
  library(xtable)
  library(usmap)

# Estimation
  library(fixest)
  library(lfe)
  library(PanelMatch)
  library(plm)
  library(estimatr)
  library(MatchIt)

# Seed
  set.seed(14618)
  
##################################################################################################################################################################################################################
# Functions
##################################################################################################################################################################################################################
# Custom Helper Functions
  source("custom_functions.R")
  
##################################################################################################################################################################################################################
# Data
##################################################################################################################################################################################################################
# Load Data
  load("vra_updated.RData") 

##################################################################################################################################################################################################################
# Generate Regional Indicators
##################################################################################################################################################################################################################
# 11 States in the former confederacy
  longdat$confed.south <- ifelse(longdat$state_name %in% c("South Carolina", "Mississippi", "Florida", "Alabama",
                                                           "Georgia", "Louisiana", "Texas", "Virginia", "Arkansas", 
                                                           "North Carolina", "Tennessee"), 1, 0)
  
  longdat$cen.south <- ifelse(longdat$state_name %in% c("Delaware", "Florida", "Georgia", "Maryland", "North Carolina",
                                                        "South Carolina", "Virginia", "West Virginia", "Alabama",
                                                        "Kentucky", "Mississippi", "Tennessee", "Arkansas", "Louisiana",
                                                        "Oklahoma", "Texas"), 1, 0)
  
  
  longdat$slave.south <- ifelse(longdat$state_name %in% c("Delaware", "Georgia", "Maryland", "South Carolina",
                                                          "Virginia", "North Carolina", "Kentucky", "Tennessee",
                                                          "Louisiana", "Mississippi", "Alabama", "Missouri", "Arkansas",
                                                          "Florida", "Texas", "West Virginia"), 1, 0)
  
##################################################################################################################################################################################################################
# Data Prep and Cleaning
##################################################################################################################################################################################################################
# Create indicator for VRA passage  
  longdat$vra <- ifelse(longdat$year > 1950, 1, 0)
  
# Clean out duplicated FIPS Codes 16087 and 51780 
  # Just drop 51780; this is Boston and South Norfolk but no elected officials outcome available for this FIPS anyway
    longdat <- longdat[which(!(longdat$fips == "51780")),]
  
  # Drop Yellowstone  
    longdat <- longdat[which(!(longdat$county_name == "Yellowstone National Park in Idaho")),]
    
  # Generate an indicator for just the original 1965 preclearance set with everyone else as control (to keep post-Allen results comparable)
    opc <- longdat[which(longdat$year == 1960), c("fips", "state_name", "county_name", "preclearance")]
    colnames(opc)[length(colnames(opc))] <- "preclearance65"
    
    longdat <- merge(longdat, opc[ , c("fips", "preclearance65")], by = "fips", all.x = T)
    
  # Generate an indicator for expansion areas
    expc <- longdat[which(longdat$year == 1970), c("fips", "preclearance")]
    colnames(expc)[length(colnames(expc))] <- "preclearance77"
    
    expansions <- merge(opc, expc, by = c("fips"), all.x = T, all.y = T)
    expansions$expansion <- as.numeric(expansions$preclearance65 != expansions$preclearance77)
    
    longdat <- merge(longdat, expansions[ , c("fips", "expansion")], by = "fips", all.x = T)
  
##################################################################################################################################################################################################################
# Main Analysis
##################################################################################################################################################################################################################
# Calculate Key Results
  # For main results I'm just going to use the 1965 original VRA covered areas and doing two versions with a 2WFE for 1957 and 1967 and directly a first difference on these two. 

# Nationwide
  # Cluster SEs and restrict everything to 57-67 for this one
    baseres <- felm(totelected ~ preclearance:vra | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960),])

  # Controls on main specification
    baseres.controls <- felm(totelected ~ preclearance:vra + I(totpop / 1000) + popdensity + pctfem + pct65plus + pctblack| factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960),])

# The South  
  # Base version with just preclearance in original coverage areas for 1965  
    baseres.south <- felm(totelected ~ preclearance:vra | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])

  # Add county control variables  
    baseres.south.controls <- felm(totelected ~ preclearance:vra + I(totpop / 1000) + popdensity + pctfem + pct65plus + pctblack| factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  
  # Now I'm going to block bootstrap those two sets of results for a comparison table on the standard errors
    baseres.south.ses.bb <- blockboot(longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,], sims = 1000, model = baseres.south, clustername = "state_name", 
                                   status = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,]$preclearance, inter = "preclearance:vra")
    
    baseres.south.controls.ses.bb <- blockbootc(longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,], sims = 1000, model = baseres.south.controls, clustername = "state_name", 
                                                status = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,]$preclearance)
    
    # SI Table 17: Compare SEs
      clusters <- as.data.frame(rbind(summary(baseres.south)$coefficients, summary(baseres.south.controls)$coefficients))
      bbs <- as.data.frame(rbind(baseres.south.ses.bb[[1]][, c(2,8)], baseres.south.controls.ses.bb[[1]][, c(2,8)]))
      comparese <- as.data.frame(cbind(clusters, bbs))
      comparese$Variable <- c("Preclearance x VRA (No Controls)", "Total Population (1,000s)", "Population per Sq. Mi.", "Pct. Female", "Pct. 65+ Years Old", "Pct. Black", "Preclearance x VRA (with Controls Above)")
      comparese <- comparese[ , c(7, 1, 2, 5, 4, 6)]
      colnames(comparese) <- c("Variable","Estimate", "SE Clustered by State", "SE Block Bootstrapped by State", "p-Value (Clustering)",  "p-Value (Block Bootstrap)")
      
      stargazer(comparese, summary = F, rownames = F, title = "Comparing Clustered and Block-Bootstrapped Standard Errors", font.size = "tiny", label = "tab:comparese", digits = 3,
                out = "si_table_17.tex",
                notes = c("Results clustered or block-bootstrapped at the state level."))
      
  # SI Table 12
      stargazer(baseres.south, baseres.south.controls, digits = 2, title = "Effect of Preclearance on County Elected Officials in the South",
                out = "si_table_12.tex", label = "tab:mainres_south",
                table.placement = "H", dep.var.labels = "County Elected Officials", 
                notes = c("Standard errors clustered by state."),
                covariate.labels = c("Total Population (1,000s)", "Pop. / Sq. Mi.", "\\% Female", "\\% 65 +", "\\% Black", "Preclearance x VRA"),
                add.lines = list(c("County and Year Fixed Effects", "\\checkmark", "\\checkmark"),
                                 c("County Demographic Controls", "", "\\checkmark"),
                                 c("Years Included", "1957, 1967", "1957, 1967")), 
                omit.stat = c("ser", "adj.rsq"),
                no.space = T, font.size = "footnotesize")
      
  # Main results figure
      baseres.south.figdat <- as.data.frame(rbind(summary(baseres.south)$coefficients["preclearance:vra",], summary(baseres.south.controls)$coefficients["preclearance:vra",]))
      baseres.south.figdat$vars <- c("Base Results", "With Controls")
      colnames(baseres.south.figdat) <- c("coef", "se", "tval", "pval", "vars")
      
      baseres.south.fig <- ggplot(data = baseres.south.figdat, aes(x = vars, y = coef)) + 
        geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
        geom_errorbar(aes(x = vars, ymin = coef - qnorm(0.95) * se , ymax = coef + qnorm(0.95) * se, width = 0.1, alpha = 1)) + 
        geom_errorbar(aes(x = vars, ymin = coef - qnorm(0.975) * se , ymax = coef + qnorm(0.975) * se, width = 0.1, alpha = 0.7)) + 
        geom_point(aes(x = vars, y = coef), size = 3) + 
        ylab("Difference-in-Differences Estimate") + xlab("") +
        geom_text(label = as.character(round(baseres.south.figdat$coef, 2)), nudge_x = 0.1) +
        ggtitle("Effect of Preclearance on County Elected Officials in the South 1957-1967") + 
        theme_bw() +
        theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
      
      ggsave(plot = baseres.south.fig, filename = "ms_figure_2.pdf", width = 7, height = 5)
      
  # Directly do Alternative Specification of this where I just first difference and run difference on preclearance dummy
    fd50 <- longdat[which(longdat$year == 1950), c("fips", "county_name", "state_name", "confed.south",  "totelected", 
                                                   "cen.south", "slave.south", "pctblack")]
    colnames(fd50)[colnames(fd50) == "pctblack"] <- "pctblack1950"
    fd60 <- longdat[which(longdat$year == 1960), c("fips", "county_name", "state_name", "preclearance", "totpop", "area", "popdensity", "pctfem", "pct65plus", "pctblack", "totelected")]
    colnames(fd60)[colnames(fd60) == "pctblack"] <- "pctblack1960"
    colnames(fd60)[colnames(fd60) == "totelected"] <- "totelected1967"
    
    fd <- merge(fd50, fd60, by = c("fips", "county_name", "state_name"), all.x = T, all.y = T)
    
    fd$elecdiff <- fd$totelected1967 - fd$totelected
    
    baseres.south.fd <- lm(elecdiff ~ preclearance, data = fd[which(fd$confed.south == 1),])
    baseres.south.fd.ses <- se(lm_robust(elecdiff ~ preclearance, data = fd[which(fd$confed.south == 1),], clusters = state_name))
    
    ptblackpop.south.fd <- lm(elecdiff ~ preclearance * pctblack1960, data = fd[which(fd$confed.south == 1),])
    ptblackpop.south.fd.ses <- se(lm_robust(elecdiff ~ preclearance * pctblack1960, data = fd[which(fd$confed.south == 1),], clusters = state_name))
    
    baseres.south.controls.fd <- lm(elecdiff ~ preclearance + I(totpop / 1000) + popdensity + pctfem + pct65plus + pctblack1960, data = fd[which(fd$confed.south == 1),])
    baseres.south.controls.fd.ses <- se(lm_robust(elecdiff ~ preclearance + I(totpop / 1000) + popdensity + pctfem + pct65plus + pctblack1960, data = fd[which(fd$confed.south == 1),], clusters = state_name))
    
    ptblackpop.south.controls.fd <- lm(elecdiff ~ preclearance * pctblack1960 +  I(totpop / 1000) + popdensity + pctfem + pct65plus, data = fd[which(fd$confed.south == 1),])
    ptblackpop.south.controls.fd.ses <- se(lm_robust(elecdiff ~ preclearance * pctblack1960 +  I(totpop / 1000) + popdensity + pctfem + pct65plus, data = fd[which(fd$confed.south == 1),], clusters = state_name))
    
    bsfd <- baseres.south.fd; bsfdse <- baseres.south.fd.ses; bpsfd <- ptblackpop.south.fd; bpsfdse <- ptblackpop.south.fd.ses; bscfd <- baseres.south.controls.fd; bscfdse <- baseres.south.controls.fd.ses
    bpscfd <- ptblackpop.south.controls.fd; bpscfdse <- ptblackpop.south.controls.fd.ses
    
    stargazer(bsfd, bscfd, bpsfd, bpscfd,
              se = list(bsfdse, bscfdse, bpsfdse, bpscfdse),
              font.size = "scriptsize",
              no.space = T,
              label = "tab:fdres",
              table.placement = "h!",
              title = "First Difference Results for Effect of Preclearance on Change in County Elected Officials",
              out =  "si_table_10.tex",
              dep.var.labels = "Change in Elected Officials 1957-67",
              notes = c("Standard errors clustered by state."),
              covariate.labels = c("Preclearance", "Total Population (1,000s)", "Pop. / Sq. Mi.", "\\% Female", "\\% 65 +", "\\% Black", "Preclearance x \\% Black", "Intercept"))
      
    # Plot a marginal effect here
      x <- seq(0, 100, 5)
      margeff.pc <- summary(ptblackpop.south.controls.fd)$coefficients["preclearance", 1] + summary(ptblackpop.south.controls.fd)$coefficients["preclearance:pctblack1960", 1] * x
      se.pc <- sqrt(vcov(ptblackpop.south.controls.fd)["preclearance", "preclearance"] + x ^ 2 * vcov(ptblackpop.south.controls.fd)["preclearance:pctblack1960", "preclearance:pctblack1960"] + 2 * x * vcov(ptblackpop.south.controls.fd)["preclearance", "preclearance:pctblack1960"])
      lb95.pc <- margeff.pc - qnorm(0.975) * se.pc
      ub95.pc <- margeff.pc + qnorm(0.975) * se.pc
      lb90.pc <- margeff.pc - qnorm(0.95) * se.pc
      ub90.pc <- margeff.pc + qnorm(0.95) * se.pc
      
      plotdat <- data.frame(cbind(x, margeff.pc, se.pc, lb90.pc, ub90.pc, lb95.pc, ub95.pc))
    
      margeff.fig <- ggplot(plotdat, aes(x = x, y = margeff.pc)) + 
        geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
        geom_errorbar(aes(x = x, ymin = lb90.pc , ymax = ub90.pc, width = 1, alpha = 1)) + 
        geom_errorbar(aes(x = x, ymin = lb95.pc, ymax = ub95.pc, width = 1, alpha = 0.7)) + 
        geom_point(aes(x = x, y = margeff.pc), size = 3) + 
        ylab("Change in County Elected Officials 1957-1967") + xlab("County % Black 1960") +
        ggtitle("Marginal Effect of Preclearance on County Elected Officials in the South 1957-67 \n by County Black Population in 1960") + 
        theme_bw() +
        theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
      
      ggsave(plot = margeff.fig, filename = "si_figure_11.pdf", width = 7, height = 5)
        
##################################################################################################################################################################################################################
# Interaction with Black Population
##################################################################################################################################################################################################################
# Merge in 1950 and 1960 black populations to fix these across counties
  pctblack50 <- longdat[longdat$year == 1950 , c("fips", "pctblack")]
  colnames(pctblack50)[2] <- "pctblack50"
  
  longdat <- left_join(longdat, pctblack50, by = "fips")
  
  pctblack60 <- longdat[longdat$year == 1960 , c("fips", "pctblack")]
  colnames(pctblack60)[2] <- "pctblack60"
  
  longdat <- left_join(longdat, pctblack60, by = "fips")
  
# Interaction with pre-treatment black population
  ptblackpop60 <- felm(totelected ~ preclearance:vra:pctblack60 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  ptblackpop60.controls <- felm(totelected ~ preclearance:vra:pctblack60 + I(totpop / 1000) + popdensity + pctfem + pct65plus + pctblack | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  
  baseres.inter <- felm(totelected ~ preclearance:vra:pctblack60 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960),])
  baseres.inter.controls <- felm(totelected ~ preclearance:vra:pctblack60 + I(totpop / 1000) + popdensity + pctfem + pct65plus + pctblack| factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960),])
  
  # All Areas Results Table      
  stargazer(baseres, baseres.controls, baseres.inter, baseres.inter.controls, digits = 2, title = "Difference-in-Differences Results Including Non-Southern Counties",
            out = "si_table_13.tex", label = "tab:mainres_all",
            table.placement = "h!", dep.var.labels = "County Elected Officials", 
            notes = c("Standard errors clustered by state."),
            covariate.labels = c("Total Population (1,000s)", "Pop. / Sq. Mi.", "\\% Female", "\\% 65 +", "\\% Black", "Preclearance x VRA", "Preclearance x VRA x 1960 Black Population"),
            add.lines = list(c("County and Year Fixed Effects", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("County Demographic Controls", "", "\\checkmark", "", "\\checkmark")), single.row = T,
            no.space = T, font.size = "footnotesize")
  
  stargazer(ptblackpop60, ptblackpop60.controls, digits = 2, title = "Additional Effect of Preclearance on Elected Officials in Counties with Larger Black Population Shares",
            out = "si_table_14.tex", label = "tab:blackpopres",
            table.placement = "H", dep.var.labels = "County Elected Officials", 
            notes = c("Standard errors clustered by state."),
            covariate.labels = c("Total Population (1,000s)", "Pop. / Sq. Mi.", "\\% Female", "\\% 65 +", "\\% Black", "Preclearance x VRA x 1960 Black Population"),
            add.lines = list(c("County and Year Fixed Effects", "\\checkmark", "\\checkmark"),
                             c("County Demographic Controls", "", "\\checkmark")), single.row = T,
            no.space = T, font.size = "footnotesize")
  
# Black population interaction results figure
  blackpop.res.south.figdat <- as.data.frame(rbind(summary(ptblackpop60)$coefficients["preclearance:vra:pctblack60",], summary(ptblackpop60.controls)$coefficients["preclearance:vra:pctblack60",]))
  blackpop.res.south.figdat$vars <- c("Base Results", "With Controls")
  colnames(blackpop.res.south.figdat) <- c("coef", "se", "tval", "pval", "vars")
  
  blackpop.res.south.fig <- ggplot(data = blackpop.res.south.figdat, aes(x = vars, y = coef)) + 
    geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
    geom_errorbar(aes(x = vars, ymin = coef - qnorm(0.95) * se , ymax = coef + qnorm(0.95) * se, width = 0.1, alpha = 1)) + 
    geom_errorbar(aes(x = vars, ymin = coef - qnorm(0.975) * se , ymax = coef + qnorm(0.975) * se, width = 0.1, alpha = 0.7)) + 
    geom_point(aes(x = vars, y = coef), size = 3) + 
    ylab("Difference-in-Differences Estimate") + xlab("") +
    geom_text(label = as.character(round(blackpop.res.south.figdat$coef, 2)), nudge_x = 0.1) +
    ggtitle("Additional Effect of Preclearance on Elected Officials in \n Counties with Larger Black Population Shares") + 
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
  
  ggsave(plot = blackpop.res.south.fig, filename = "ms_figure_3.pdf", width = 7, height = 5)  
  
  
# Merge in Matthews and Prothro Data and do this with registered voters
# Load Matthews and Prothro Data on the South (http://www.icpsr.umich.edu/icpsrweb/ICPSR/studies/7255)
# Original data has states but no concordance for county codes. County names merged in assuming counties were just
# listed in alphabetical order
  mp <- read.csv("matpro.csv", header = T, stringsAsFactors = F)
  
# Harmonize fips  
  mp$fips <- ifelse(nchar(mp$fips) == 4, paste("0", mp$fips, sep = ""), mp$fips)
  
# Just grab columns I care about to avoid blowing up panel and repeat columns  
  colnames(mp)[colnames(mp) == "state.abb"] <- "mp.state.abb"
  colnames(mp)[colnames(mp) == "county.name"] <- "mp.county.name"
  mp <- mp[ , c("fips", "mp.state.abb", "mp.county.name", "bvap60", "wvap60","pbvprg60", "pbpreg64", "pbpreg67", "pbvprg68", "pwvprg60", "pwpreg64", "pwpreg67", "pwvprg68", "fedex66", "fedexs69")]
  
# Merge data
  longdat <- merge(longdat, mp, by = "fips", all.x = T)
 
# Just VAP?
  bvap <- felm(totelected ~ preclearance:vra:bvap60 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  
# VAP As Share of Total Electorate
  longdat$bvapshare <- longdat$bvap60 / (longdat$bvap60 + longdat$wvap60) * 100
  
  bvapshare <- felm(totelected ~ preclearance:vra:bvapshare | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  
# What about counties where black registration already growing (and scaring incumbents)? Nothing here
  longdat$blackreggrowth <- longdat$pbpreg67 - longdat$pbvprg60  
  
  blackreggrowth <- felm(totelected ~ preclearance:vra:blackreggrowth | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  
# Examiners?
  examiners <- felm(totelected ~ vra:fedex66 | factor(fips) + factor(year) | 0 | state_name, data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  examiners69 <- felm(totelected ~ vra:fedexs69 | factor(fips) + factor(year) | 0 | state_name, data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  
  stargazer(bvap, bvapshare, blackreggrowth, digits = 2, title = "Additional Effect of Preclearance on Elected Officials in Counties with Larger Black Voting Age or Registered Populations", 
            out = "si_table_7.tex", label = "tab:bvap_res",
            table.placement = "h!", dep.var.labels = "County Elected Officials", 
            notes = c("Standard errors clustered by state."),
            covariate.labels = c("Preclearance x VRA x 1960 BVAP", "Preclearance x VRA x Black Share of 1960 VAP", "Preclearance x VRA x Growth in \\% of BVAP Registered 1960-67"),
            single.row = T,
            no.space = T, font.size = "footnotesize")
  
  
  stargazer(examiners, digits = 2, title = "Additional Effect of Preclearance on Elected Officials in Counties with Federal Examiners", 
            out = "si_table_8.tex", label = "tab:examiners",
            table.placement = "h!", dep.var.labels = "County Elected Officials", 
            notes = c("Standard errors clustered by state."),
            covariate.labels = c("Preclearance x VRA x Federal Examiners"),
            single.row = T,
            no.space = T, font.size = "footnotesize")
  
##################################################################################################################################################################################################################
# Interaction with Democratic Vote Share
##################################################################################################################################################################################################################
# All replication files for Democratic vote shares rely on Dave Leip's Election Data. This data is proprietary and cannot be shared. I'm
# leaving the code here commented out so a user can replicate this assuming they have institutional access to the data or purchase it but
# I cannot provide the extracts I purchased in compliance with their use policy. 
  
  #votes32 <- read_csv("County_Presidential_Election_Data_1932_0_0_2.csv", skip = 1)
  #votes36 <- read_csv("County_Presidential_Election_Data_1936_0_0_2.csv", skip = 1)
  #votes40 <- read_csv("County_Presidential_Election_Data_1940_0_0_2.csv", skip = 1)
  #votes44 <- read_csv("County_Presidential_Election_Data_1944_0_0_2.csv", skip = 1)
  #votes48 <- read_csv("County_Presidential_Election_Data_1948_0_0_2.csv", skip = 1)
  #votes52 <- read_csv("County_Presidential_Election_Data_1952_0_0_2.csv", skip = 1)
  #votes56 <- read_csv("County_Presidential_Election_Data_1956_0_0_2.csv", skip = 1)
  #votes60 <- read_csv("County_Presidential_Election_Data_1960_0_0_2.csv", skip = 1)
  #votes64 <- read_csv("County_Presidential_Election_Data_1964_0_0_2.csv", skip = 1)
  
  #votes32$fips <- ifelse(nchar(votes32$fips) == 4, paste(0, votes32$fips, sep = ""), votes32$fips)
  #votes32$dshare2p32 <- votes32$vote1 / (votes32$vote1 + votes32$vote2) * 100
  
  #votes36$fips <- ifelse(nchar(votes36$fips) == 4, paste(0, votes36$fips, sep = ""), votes36$fips)
  #votes36$dshare2p36 <- votes36$vote1 / (votes36$vote1 + votes36$vote2) * 100
  
  #votes40$fips <- ifelse(nchar(votes40$fips) == 4, paste(0, votes40$fips, sep = ""), votes40$fips)
  #votes40$dshare2p40 <- votes40$vote1 / (votes40$vote1 + votes40$vote2) * 100
  
  #votes44$fips <- ifelse(nchar(votes44$fips) == 4, paste(0, votes44$fips, sep = ""), votes44$fips)
  #votes44$dshare2p44 <- votes44$vote1 / (votes44$vote1 + votes44$vote2) * 100
  
  #votes48$fips <- ifelse(nchar(votes48$fips) == 4, paste(0, votes48$fips, sep = ""), votes48$fips)
  #votes48$dshare2p48 <- votes48$vote1 / (votes48$vote1 + votes48$vote2) * 100
  
  #votes52$fips <- ifelse(nchar(votes52$fips) == 4, paste(0, votes52$fips, sep = ""), votes52$fips)
  #votes52$dshare2p52 <- votes52$vote1 / (votes52$vote1 + votes52$vote2) * 100
  
  #votes56$fips <- ifelse(nchar(votes56$fips) == 4, paste(0, votes56$fips, sep = ""), votes56$fips)
  #votes56$dshare2p56 <- votes56$vote1 / (votes56$vote1 + votes56$vote2) * 100
  
  #votes60$fips <- ifelse(nchar(votes60$fips) == 4, paste(0, votes60$fips, sep = ""), votes60$fips)
  #votes60$dshare2p60 <- votes60$vote1 / (votes60$vote1 + votes60$vote2) * 100
  
  #votes64$fips <- ifelse(nchar(votes64$fips) == 4, paste(0, votes64$fips, sep = ""), votes64$fips)
  #votes64$dshare2p64 <- votes64$vote1 / (votes64$vote1 + votes64$vote2) * 100
  
  #votes <- merge(votes32[, c("fips", "name", "type", "dshare2p32")], votes36[ , c("fips", "dshare2p36")], by = "fips", all.x = T)
  #votes <- merge(votes, votes40[ , c("fips", "dshare2p40")], by = "fips", all.x  = T)
  #votes <- merge(votes, votes44[ , c("fips", "dshare2p44")], by = "fips", all.x  = T)
  #votes <- merge(votes, votes48[ , c("fips", "dshare2p48")], by = "fips", all.x  = T)
  #votes <- merge(votes, votes52[ , c("fips", "dshare2p52")], by = "fips", all.x  = T)
  #votes <- merge(votes, votes56[ , c("fips", "dshare2p56")], by = "fips", all.x  = T)
  #votes <- merge(votes, votes60[ , c("fips", "dshare2p60")], by = "fips", all.x  = T)
  #votes <- merge(votes, votes64[ , c("fips", "dshare2p64")], by = "fips", all.x  = T)
  
  #longdat <- merge(longdat, votes[ , -c(2,3)], by = "fips", all.x = T)
  
# Regression Results
  #most.dem.1932 <- felm(totelected ~ preclearance:vra:dshare2p32 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  #most.dem.1932.controls <- felm(totelected ~ preclearance:vra:dshare2p32 + + I(totpop / 1000) + popdensity + pctfem + pct65plus + pctblack | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  #most.dem.1936 <- felm(totelected ~ preclearance:vra:dshare2p36 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  #most.dem.1940 <- felm(totelected ~ preclearance:vra:dshare2p40 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  #most.dem.1944 <- felm(totelected ~ preclearance:vra:dshare2p44 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  #most.dem.1948 <- felm(totelected ~ preclearance:vra:dshare2p48 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  #most.dem.1952 <- felm(totelected ~ preclearance:vra:dshare2p52 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  #most.dem.1956 <- felm(totelected ~ preclearance:vra:dshare2p56 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  #most.dem.1960 <- felm(totelected ~ preclearance:vra:dshare2p60 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  #most.dem.1964 <- felm(totelected ~ preclearance:vra:dshare2p64 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$confed.south == 1,])
  
  #demshare.res.south.figdat <- as.data.frame(rbind(summary(most.dem.1932)$coefficients["preclearance:vra:dshare2p32",], summary(most.dem.1932.controls)$coefficients["preclearance:vra:dshare2p32",]))
  #demshare.res.south.figdat$vars <- c("Base Results", "With Controls")
  #colnames(demshare.res.south.figdat) <- c("coef", "se", "tval", "pval", "vars")
  
  #demshare.res.south.fig <- ggplot(data = demshare.res.south.figdat, aes(x = vars, y = coef)) + 
  #  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  #  geom_errorbar(aes(x = vars, ymin = coef - qnorm(0.95) * se , ymax = coef + qnorm(0.95) * se, width = 0.1, alpha = 1)) + 
  #  geom_errorbar(aes(x = vars, ymin = coef - qnorm(0.975) * se , ymax = coef + qnorm(0.975) * se, width = 0.1, alpha = 0.7)) + 
  #  geom_point(aes(x = vars, y = coef), size = 3) + 
  #  ylab("Difference-in-Differences Estimate") + xlab("") +
  #  geom_text(label = as.character(round(demshare.res.south.figdat$coef, 2)), nudge_x = 0.1) +
  #  ggtitle("Additional Effect of Preclearance on Elected Officials in \n Counties with Larger Democratic Vote Shares") + 
  #  theme_bw() +
  #  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
  
#  ggsave(plot = demshare.res.south.fig, filename = "party_res.pdf", width = 7, height = 5)  
  
#  stargazer(most.dem.1932, most.dem.1932.controls, digits = 2, 
#            title = "Additional Effect of Preclearance on Elected Officials in Counties with Larger Democratic Vote Shares",
#            out = "mostdem.tex", label = "tab:mostdem",
#            table.placement = "h!", dep.var.labels = "County Elected Officials", 
#            notes = c("Standard errors clustered by state."),
#            covariate.labels = c("Total Population (1,000s)", "Pop. / Sq. Mi.", "\\% Female", "\\% 65 +", "\\% Black", "Preclearance x Allen x \\% Dem. 1932"),
#            add.lines = list(c("County and Year Fixed Effects", "\\checkmark", "\\checkmark"),
#                             c("County Demographic Controls", "", "\\checkmark")), single.row = T,
#            no.space = T, font.size = "footnotesize") 
  
# Plot Across Years to Show Choice of Year Doesn't Matter
#  demshare.res.south.multiyear.figdat <- as.data.frame(rbind(summary(most.dem.1932)$coefficients["preclearance:vra:dshare2p32",], 
#                                                   summary(most.dem.1936)$coefficients["preclearance:vra:dshare2p36",],
#                                                   summary(most.dem.1940)$coefficients["preclearance:vra:dshare2p40",],
#                                                   summary(most.dem.1944)$coefficients["preclearance:vra:dshare2p44",],
#                                                   summary(most.dem.1948)$coefficients["preclearance:vra:dshare2p48",],
#                                                   summary(most.dem.1952)$coefficients["preclearance:vra:dshare2p52",],
#                                                   summary(most.dem.1956)$coefficients["preclearance:vra:dshare2p56",],
#                                                   summary(most.dem.1960)$coefficients["preclearance:vra:dshare2p60",],
#                                                   summary(most.dem.1964)$coefficients["preclearance:vra:dshare2p64",]))
#  demshare.res.south.multiyear.figdat$vars <- c("1932", "1936", "1940", "1944", "1948", "1952", "1956", "1960", "1964")
#  colnames(demshare.res.south.multiyear.figdat) <- c("coef", "se", "tval", "pval", "vars")
  
#  demshare.res.south.multiyear.fig <- ggplot(data = demshare.res.south.multiyear.figdat, aes(x = vars, y = coef)) + 
#    geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
#    geom_errorbar(aes(x = factor(vars), ymin = coef - qnorm(0.95) * se , ymax = coef + qnorm(0.95) * se, width = 0.1, alpha = 1)) + 
#    geom_errorbar(aes(x = factor(vars), ymin = coef - qnorm(0.975) * se , ymax = coef + qnorm(0.975) * se, width = 0.1, alpha = 0.7)) + 
#    geom_point(aes(x = factor(vars), y = coef), size = 3) + 
#    ylab("Difference-in-Differences Estimate") + xlab("") + 
#    geom_text(label = as.character(round(demshare.res.south.multiyear.figdat$coef, 2)), nudge_x = 0.3, angle = 45) +
#    ggtitle("Additional Effect of Preclearance on Elected Officials in \n Counties with Larger Democratic Vote Shares") + 
#    theme_bw() +
#    theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
  
#  ggsave(plot = demshare.res.south.multiyear.fig, filename = "party_res_years.pdf", width = 7, height = 5)  
  
  
##################################################################################################################################################################################################################
# North Carolina
##################################################################################################################################################################################################################
  nc <- longdat[which(longdat$state_name == "North Carolina"),]

  nc.main <- felm(totelected ~ preclearance:vra | factor(fips) + factor(year), data = nc[nc$year %in% c(1950, 1960),])

  nc.inter <- felm(totelected ~ preclearance:vra:pctblack60 | factor(fips) + factor(year),  data = nc[nc$year %in%  c(1950, 1960),])

  #nc.party <- felm(totelected ~ preclearance:vra:dshare2p32 | factor(fips) + factor(year),  data = nc[nc$year %in%  c(1950, 1960),])

  nc.bvap <- felm(totelected ~ preclearance:vra:bvap60 | factor(fips) + factor(year),  data = nc[nc$year %in%  c(1950, 1960),])

# Version with party  
  #stargazer(nc.main, nc.inter, nc.bvap, nc.party, digits = 2, title = "North Carolina Difference-in-Differences Results", 
  #        out = "nc.tex", label = "tab:nc",
  #        table.placement = "H", dep.var.labels = "County Elected Officials", covariate.labels = c("Covered x VRA", "Covered x VRA x 1960 Black Pop.", "Covered x VRA x 1960 Black VAP", "Covered x VRA x 1932 Dem. Vote Share"),
  #        no.space = T, font.size = "scriptsize")
  
# Version without party  
  stargazer(nc.main, nc.inter, nc.bvap, digits = 2, title = "North Carolina Difference-in-Differences Results", 
          out = "ms_table_4.tex", label = "tab:nc",
          table.placement = "H", dep.var.labels = "County Elected Officials", covariate.labels = c("Covered x VRA", "Covered x VRA x 1960 Black Pop.", "Covered x VRA x 1960 Black VAP", "Covered x VRA x 1932 Dem. Vote Share"),
          no.space = T, font.size = "scriptsize")
  
##################################################################################################################################################################################################################
# After the Allen Decision
##################################################################################################################################################################################################################
# Post - Allen
  longdat$allen <- ifelse(longdat$year > 1960, 1, 0)
  
# Include all the post-Allen years but only do this for the set of counties that were preclearance originally as treatment. Set the remaining ones to control
  # Do This with Just 1965 Treated
  postallen.65set.base <- felm(totelected ~ preclearance65:allen | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$confed.south == 1 & longdat$year %in% c(1960, 1970, 1980, 1990) & longdat$expansion == 0,])
  
  postallen.65set.inter <- felm(totelected ~ preclearance65:allen:pctblack60 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$confed.south == 1 & longdat$year %in% c(1960, 1970, 1980, 1990) & longdat$expansion == 0,])
  
  postallen.65set.controls <- felm(totelected ~ preclearance65:allen + I(totpop / 1000) + popdensity + pctfem + pct65plus + pctblack| factor(fips) + factor(year)| 0 | state_name,  data = longdat[longdat$year %in% c(1960, 1970, 1980, 1990) & longdat$confed.south == 1 & longdat$expansion == 0,])
  
  postallen.65set.inter.controls <- felm(totelected ~ preclearance65:allen:pctblack60 + I(totpop / 1000) + popdensity + pctfem + pct65plus + pctblack | factor(fips) + factor(year)| 0 | state_name,  data = longdat[longdat$confed.south == 1 & longdat$year %in% c(1960, 1970, 1980, 1990) & longdat$expansion == 0,])
  
  postallen.65set.figdat <- as.data.frame(rbind(summary(postallen.65set.base)$coefficients["preclearance65:allen",], 
                                                summary(postallen.65set.controls)$coefficients["preclearance65:allen",]))
  postallen.65set.figdat$vars <- c("Base", "Controls")
  colnames(postallen.65set.figdat) <- c("coef", "se", "tval", "pval", "vars")
  
  postallen.65set.fig <- ggplot(data = postallen.65set.figdat, aes(x = vars, y = coef)) + 
    geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
    geom_errorbar(aes(x = factor(vars), ymin = coef - qnorm(0.95) * se , ymax = coef + qnorm(0.95) * se, width = 0.1, alpha = 1)) + 
    geom_errorbar(aes(x = factor(vars), ymin = coef - qnorm(0.975) * se , ymax = coef + qnorm(0.975) * se, width = 0.1, alpha = 0.7)) + 
    geom_point(aes(x = factor(vars), y = coef), size = 3) + 
    ylab("Difference-in-Differences Estimate") + xlab("") + 
    geom_text(label = as.character(round(postallen.65set.figdat$coef, 2)), nudge_x = 0.1) +
    ggtitle("Elected Officials in Covered Counties After Allen (1969)") + 
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
  
  ggsave(plot = postallen.65set.fig, filename = "ms_figure_5.pdf", width = 7, height = 5) 
  
  postallen.65set.inter.figdat <- as.data.frame(rbind(summary(postallen.65set.inter)$coefficients["preclearance65:allen:pctblack60",], 
                                                summary(postallen.65set.inter.controls)$coefficients["preclearance65:allen:pctblack60",]))
  postallen.65set.inter.figdat$vars <- c("Base", "Controls")
  colnames(postallen.65set.inter.figdat) <- c("coef", "se", "tval", "pval", "vars")
  
  postallen.65set.inter.fig <- ggplot(data = postallen.65set.inter.figdat, aes(x = vars, y = coef)) + 
    geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
    geom_errorbar(aes(x = factor(vars), ymin = coef - qnorm(0.95) * se , ymax = coef + qnorm(0.95) * se, width = 0.1, alpha = 1)) + 
    geom_errorbar(aes(x = factor(vars), ymin = coef - qnorm(0.975) * se , ymax = coef + qnorm(0.975) * se, width = 0.1, alpha = 0.7)) + 
    geom_point(aes(x = factor(vars), y = coef), size = 3) + 
    ylab("Difference-in-Differences Estimate") + xlab("") + 
    geom_text(label = as.character(round(postallen.65set.inter.figdat$coef, 2)), nudge_x = 0.1) +
    ggtitle("Elected Officials in Covered Counties with Larger Black \n Populations After Allen (1969)") + 
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
  
  ggsave(plot = postallen.65set.inter.fig, filename = "ms_figure_6.pdf", width = 7, height = 5) 
  
  pab <- postallen.65set.base; pac <- postallen.65set.controls; paib <- postallen.65set.inter; paic <- postallen.65set.inter.controls
  
  stargazer(pab, pac, paib, paic, digits = 2, 
            title = "Effect of Preclearance on County Elected Officials After Allen v. State Board of Elections (1969)",
            out = "si_table_16.tex", label = "tab:allen_65set",
            table.placement = "h!", dep.var.labels = "County Elected Officials", 
            notes = c("Results including only never-covered counties and counties originally covered in 1965.",  "Standard errors clustered by state."),
            covariate.labels = c("Total Population (1,000s)", "Pop. / Sq. Mi.", "\\% Female", "\\% 65 +", "\\% Black", "Preclearance x Allen", "Preclearance x Allen x \\% Black"),
            add.lines = list(c("County and Year Fixed Effects", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("County Demographic Controls", "", "\\checkmark", "", "\\checkmark")), single.row = T,
            no.space = T, font.size = "footnotesize")
  
# For robustness just include everything with true treatment status for the year
  postallen.fullset.base <- felm(totelected ~ preclearance:allen | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$confed.south == 1 & longdat$year %in% c(1960, 1970, 1980, 1990),])
  
  postallen.fullset.inter <- felm(totelected ~ preclearance:allen:pctblack60 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$confed.south == 1 & longdat$year %in% c(1960, 1970, 1980, 1990),])
  
  postallen.fullset.controls <- felm(totelected ~ preclearance:allen + I(totpop / 1000) + popdensity + pctfem + pct65plus + pctblack| factor(fips) + factor(year)| 0 | state_name,  data = longdat[longdat$year %in% c(1960, 1970, 1980, 1990) & longdat$confed.south == 1,])
  
  postallen.fullset.inter.controls <- felm(totelected ~ preclearance:allen:pctblack60 + I(totpop / 1000) + popdensity + pctfem + pct65plus + pctblack | factor(fips) + factor(year)| 0 | state_name,  data = longdat[longdat$confed.south == 1 & longdat$year %in% c(1960, 1970, 1980, 1990),])
  
  pafb <- postallen.fullset.base; pafc <- postallen.fullset.controls; paifb <- postallen.fullset.inter; paifc <- postallen.fullset.inter.controls
  
  stargazer(pafb, pafc, paifb, paifc, digits = 2, 
            title = "Effect of Preclearance on County Elected Officials After Allen v. State Board of Elections (1969)",
            out = "si_table_11.tex", label = "tab:allen_fullset",
            table.placement = "h!", dep.var.labels = "County Elected Officials", 
            notes = c("Standard errors clustered by state."),
            covariate.labels = c("Total Population (1,000s)", "Pop. / Sq. Mi.", "\\% Female", "\\% 65 +", "\\% Black", "Preclearance x Allen", "Preclearance x Allen x \\% Black"),
            add.lines = list(c("County and Year Fixed Effects", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("County Demographic Controls", "", "\\checkmark", "", "\\checkmark")), single.row = T,
            no.space = T, font.size = "footnotesize")

##################################################################################################################################################################################################################
# Summary Statistics
##################################################################################################################################################################################################################
# Officials and County Coverage by Year
  out <- data.table(longdat)
  out2 <- out[is.na(out$totelected) == 0 , list(counties = length(fips),
                                                elected = sum(totelected, na.rm = T)), by = list(year)]
  
  setkey(out2, year)
  out2$year <- c("1957", "1967", "1977", "1987", "1992")
  colnames(out2) <- c("Year", "Counties", "Elected Officials")
  
  stargazer(out2, summary = F, font.size = "scriptsize", rownames = F, digits = 0, no.space = T,
            title = "Nationwide Counties and Elected Officials Covered by Year",
            out = "ms_table_2.tex",
            label = "tab:cogcoverage", table.placement = "H")  
  
# Now We Can Look at Summary Statistics by State
  state.diffs <- fd %>% filter(!is.na(elecdiff) & confed.south == 1) %>% 
    group_by(state_name, preclearance) %>%
    summarize(countiesindata = length(elecdiff),
              coveredcounties = length(elecdiff[preclearance == 1]),
              electeds = round(median(totelected, na.rm = T), 0),
              avgchange = round(mean(elecdiff), 2),
              losscount = round(sum(elecdiff < 0) / length(elecdiff), 2),
              avgpctblack = round(mean(pctblack1950, na.rm = T), 2),
              preclearance = round(ifelse(sum(preclearance) > 0, 1,0), 0)) %>% 
    arrange(preclearance)
  
  
  colnames(state.diffs) <- c("State", "Preclearance", "Total Counties", "Covered Counties", "Median Elected Officials 1950","Avg. Change 1957-67", "% Counties w/ Negative Avg. Change", 
                             "Avg. % Black 1950")
  
  stargazer(state.diffs, summary = F, font.size = "tiny", rownames = F, digits = 2,
            column.labels =  c("State", "Preclearance", "Total Counties", "Covered Counties", "Median Elected Officials 1950","Avg. Diff. in Officials 1957-67", "% of Counties where Avg. Diff. in Officials 1957-67 < 0", 
                               "Avg. County % Black 1950"),
            title = "VRA Coverage, Elected Officials, and Demographic Composition in Southern Counties",
            out = "ms_table_3.tex",
            label = "tab:statesumm", table.placement = "h!",
            notes = c("Note: Counties refers to the number of counties in each state for which data on elected officials is available. This does not include consolidated city-county governments or",
                      "independent cities which perform many county government functions (as in Virginia). Avg. Change 1957-67 refers to average county-level change in total elected officials",
                      "1957-67. \\% Counties w/ Avg. Change refers to the proportion of counties within each state that lost elected officials between 1957 and 1967."),
            column.sep.width = "-7pt")
  
  
##################################################################################################################################################################################################################
# Schools
##################################################################################################################################################################################################################
# complied manually from Census of Governments qualitative notes on schools by state  
  out <- data.frame(State = state.name, Counted = rep("NA", length(state.name)), Offices = rep("NA", length(state.name)))
  out$Counted <- c("N", "NA", "Y", "Y", "Y", "Y", "N", "N", "N", "N", "NA", "N", "Y", "N", "Y", "Y", "N", "N", "N", "Y", "N", "N", "Y", "Y", "Y", "Y", "Y",
                   "N", "N", "N", "N", "N", "N", "Y", "Y", "Y", "Y", "N", "NA", "Y", "Y", "Y", "Y", "N", "N", "N", "Y", "N", "Y", "Y")
  out$Offices <- c("NA", "NA", "Superintendent", "BOE", "BOE, Superintendent", "Superintendent", "NA", "NA", "NA", "NA", "NA","NA", "BOE, Land Assessor for Schools, Superintendent",
                   "NA", "BOE", "BOE, Superintendent", "NA", "NA", "NA", "BOE", "NA", "NA", "Superintendent", "BOE, Superintendent", "Superintendent", "Superintendent", "Superintendent",
                   "NA", "NA", "NA", "NA", "NA", "NA", "Superintendent", "BOE", "Superintendent", "BOE", "NA", "NA", "BOE, Superintendent", "BOE, Superintendent", "BOE, Superintendent",
                   "BOE, Superintendent", "NA", "NA", "NA", "BOE, Superintendent", "NA", "Superintendent", "Superintendent")
  out$Notes <- c("NA", "Not a state until 1959", rep("NA", 8), "Not a state until 1959",rep("NA", 8), "Montgomery County Only", rep("NA", 16), "Rural Schools Only", rep("NA", 13))
  out$change67 <- c("N", "NA", rep("N", 4), "County governments abolished in 1960", "N","N", "N", "NA", rep("N", 4), "BOE removed in 1967", rep("N", 7), "BOE removed in 1967", "BOE added in 1967", rep("N", 11), "Only counties without independent school districts", rep("N", 11), "Superintendent removed in 1967", "N")
  colnames(out) <- c("State", "Counted", "Offices", "Notes", "Changes to County Offices 1957-67")
  
  stargazer(out, summary = F, font.size = "tiny", rownames = F, digits = 0, notes.align = "l",
            title = "Education-Related Elected Offices Included in County Governments by State",
            out = "si_table_6.tex",
            label = "tab:stateedsumm", table.placement = "H",
            notes = c("Note: BOE is shorthand used to indicate ``Board of Education,'' or school governing body, even in states where these go by by different names. A Y in the Counted column indicates",
                      "that the Census counts any school-related offices as a part of county governments. The notes column identifies states where BOEs are elected but only in specific counties. All data",
                      "summarized here was manually transcribed from the  1957 Census of Governments."))
  
##################################################################################################################################################################################################################
# Parallel Trends in Elected Officials
##################################################################################################################################################################################################################
# Trends in Elected Officials  
  # Generate preclearance label that identifies the counties that would be covered but are not in 1957
    longdat$preclearancelabel <- longdat$preclearance
  
    out <- data.table(longdat)
    out <- out[is.na(preclearance) == 0 & year == 1960, list(pc = mean(preclearance, na.rm = T), state_name = state_name), by = list(fips)]  
  
    longdat <- merge(longdat, out[, c("fips", "pc")], by = "fips", all.x = T)
  
    longdat$preclearancelabel[longdat$year == 1950] <- ifelse(longdat$pc[longdat$year == 1950] > 0, 1, 0)
  
    out <- data.table(longdat)
    out2 <- out[is.na(preclearance) == 0 & confed.south == 1, list(totelected = mean(totelected, na.rm = T)), by = list(preclearancelabel, year)]
    out2$yrlabel <- ifelse(out2$year == 1950, 1957,
                           ifelse(out2$year == 1960, 1967,
                                  ifelse(out2$year == 1970, 1977,
                                         ifelse(out2$year == 1980, 1987, 1992))))
    
    electtrends.south <- ggplot(out2, aes(x = factor(yrlabel), y = totelected, group = factor(preclearancelabel))) + 
      geom_line(data = out2, aes(linetype = factor(preclearancelabel))) + scale_x_discrete(breaks = c(1957, 1967, 1977, 1987, 1992)) +
      geom_vline(xintercept = 1.8, color = "red") + geom_vline(xintercept = 2.2, color = "red") + ggtitle("County Elected Officials") + xlab("") + 
      annotate("text", x = 1.2, y = 30, label = "Voting Rights Act \n (1965)", color = "red") +
      annotate("text", x = 3, y = 25.25, label = "Allen v. State Board of \n Elections (1969)", color = "red") +
      ylab("Average # of Elected Officials per County") + labs(linetype = "Preclearance") + theme_bw() + 
      theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
    
    ggsave(plot = electtrends.south, filename = "ms_figure_1_left.pdf", width = 7, height = 5)
    
  
##################################################################################################################################################################################################################
# Parallel Trends with General Control and Police Employees
##################################################################################################################################################################################################################
# Parallel Trends
# 1957
  # Grab Public Employment from 1957
  # Financial administration rolled into general control in 1957 from the documentation
  pubemp57 <- read.csv("publicemps.csv", header = T, stringsAsFactors = F)
  pubemp57$state <- toupper(pubemp57$state)
  pubemp57 <- pubemp57[ , c("state", "county", "general.control", "police")]
  colnames(pubemp57) <- c("state", "county", "govemp", "police")
  pubemp57$year <- 1957

# 1962
# Grab Employees in 1962    
  p.raw <- readLines("00017-0001-Data.txt")  

# Format Separate Fields
  pubemp62 <- data.frame(type = substr(p.raw, 1, 1), 
                         icpstate = substr(p.raw, 5, 6),
                         county = trimws(substr(p.raw, 9, 25), "both"),
                         totalemp = as.numeric(substr(p.raw, 85, 90)),
                         fte = as.numeric(substr(p.raw, 91, 96)),
                         gencont = as.numeric(substr(p.raw, 176, 180)),
                         finadmin = as.numeric(substr(p.raw, 171, 175)),
                         police = as.numeric(substr(p.raw, 161, 165)))
  
  icpstates <- unique(pubemp62[pubemp62$type == "S",])
  pubemp62 <- merge(pubemp62, icpstates[, c("icpstate", "county")], by = "icpstate", all.x = T)
  pubemp62 <- pubemp62[which(pubemp62$type == "C"),]
  
  pubemp62$govemp <- pubemp62$finadmin + pubemp62$gencont
  
  pubemp62 <- pubemp62[ , c("county.y", "county.x", "govemp", "police")]
  
  colnames(pubemp62) <- c("state", "county", "govemp", "police")
  
  pubemp62$year <- 1962

# 1967
# The readlines txt didn't match the pdf here, so using the ICPSR created rda that does. Readlines was also missing police
  load("1967_cog.rda")
  p.raw67 <- da00017.0002; rm(da00017.0002)
  
  pubemp67 <- p.raw67[, c("TYPE", "YEAR", "STATE", "LOCATION", "FT_EMPLOY_POLICE", "FT_EMPLOY_FINADMIN", "FT_EMPLOY_GENCNTRL")]
  pubemp67$TYPE <- as.character(pubemp67$TYPE)
  pubemp67 <- pubemp67[which(pubemp67$TYPE == "(1) County"),]
  
  pubemp67$STATE <- as.character(pubemp67$STATE)
  pubemp67$STATE <- toupper(splitit(pubemp67$STATE, ") ", 2))
  pubemp67$LOCATION <- trimws(as.character(pubemp67$LOCATION), which = "both")
  
  pubemp67$govemp <- pubemp67$FT_EMPLOY_GENCNTRL + pubemp67$FT_EMPLOY_FINADMIN
  
  pubemp67 <- pubemp67[ , c("STATE", "LOCATION", "govemp", "FT_EMPLOY_POLICE", "YEAR")]
  
  colnames(pubemp67) <- c("state", "county", "govemp", "police", "year")

# 1977  
  p.raw77 <- readLines("08117-0004-Data.txt")  
  pubemp77 <- data.frame(statecode = substr(p.raw77, 1, 2),
                         type = trimws(substr(p.raw77, 3, 3), which = "both"),
                         name = substr(p.raw77, 10, 45),
                         finadminfull = as.numeric(substr(p.raw77, 412, 423)),
                         gencontfull = as.numeric(substr(p.raw77, 460, 471)),
                         police = as.numeric(substr(p.raw77, 604, 615)))
  
  icpstates <- unique(pubemp77[pubemp77$type == "0",])
  icpstates$name <- splitit(splitit(icpstates$name, "STATE OF ", 2), "\\.", 1)
  pubemp77 <- merge(pubemp77, icpstates[, c("statecode", "name")], by = "statecode", all.x = T)
  pubemp77 <- pubemp77[which(pubemp77$type == "1"),]
  
  pubemp77$name.x <- splitit(pubemp77$name.x, " COUNTY" , 1)
  pubemp77$name.x <- splitit(pubemp77$name.x, " BOROUGH" , 1)
  pubemp77$name.x <- splitit(pubemp77$name.x, " PARISH" , 1)
  pubemp77$govemp <- pubemp77$gencontfull + pubemp77$finadminfull
  
  pubemp77 <- pubemp77[ , c("name.y", "name.x", "govemp", "police")]
  
  colnames(pubemp77) <- c("state", "county", "govemp", "police")
  
  pubemp77$year <- 1977

# 1987  
  pubemp87 <- read_dta("1987data.dta")
  pubemp87 <- pubemp87[ , c("V1A", "V1B", "V5", "V36", "V42", "V54", "V60")]
  colnames(pubemp87) <- c("state", "type", "county", "finadminftee", "gengovftee", "policearrest", "policeother")
  pubemp87$county <- trimws(splitit(pubemp87$county, "COUNTY", 1), which = "both")
  pubemp87$county <- trimws(splitit(pubemp87$county, "BOROUGH", 1), which = "both")
  pubemp87$county <- trimws(splitit(pubemp87$county, "PARISH", 1), which = "both")
  pubemp87 <- pubemp87[which(pubemp87$type == 1),]
  
  icpstates <- data.frame(attributes(pubemp87$state)$labels)
  icpstates$state <- rownames(icpstates)
  colnames(icpstates) <- c("state", "state.name")
  
  pubemp87 <- merge(pubemp87, icpstates, by = "state", all.x = T)
  
  pubemp87$govemp <- pubemp87$finadminftee + pubemp87$gengovftee
  
  pubemp87$police <- pubemp87$policearrest + pubemp87$policeother
  
  pubemp87 <- pubemp87[ , c("state.name", "county", "govemp", "police")]
  
  colnames(pubemp87) <- c("state", "county", "govemp", "police")
  
  pubemp87$year <- 1987

#Combine
  pubemp.full <- data.frame(rbind(pubemp57, pubemp62, pubemp67, pubemp77, pubemp87))
  
  pubemp.full$pc <- ifelse(pubemp.full$state %in% toupper(c("Alabama", "Georgia", "Louisiana", "Mississippi", "South Carolina", "Virginia")), 1, 0)
  pubemp.full$pc <- ifelse(pubemp.full$state == "NORTH CAROLINA" & pubemp.full$county %in% c("ANSON", "BEAUFORT", "BERTIE", "BLADEN", "CAMDEN", "CASWELL", "CHOWAN", "CLEVELAND", "CRAVEN", "CUMBERLAND",
                                                                                             "EDGECOMBE", "FRANKLIN", "GASTON", "GATES", "GRANVILLE", "GREENE", "GUILFORD", "HALIFAX", "HARNETT", "HERTFORD",
                                                                                             "HOKE", "LEE", "LENOIR", "MARTIN", "NASH", "NORTHAMPTON", "ONSLOW", "PASQUOTANK", "PERQUIMANS", "PERSON",
                                                                                             "PITT", "ROBESON", "ROCKINGHAM", "SCOTLAND", "UNION", "VANCE", "WASHINGTON", "WAYNE", "WILSON"), 1, pubemp.full$pc)
  
  pubemp.full$confed.south <- ifelse(pubemp.full$state %in% toupper(c("South Carolina", "Mississippi", "Florida", "Alabama",
                                                                      "Georgia", "Louisiana", "Texas", "Virginia", "Arkansas", 
                                                                      "North Carolina", "Tennessee")), 1, 0)
  
  out2 <- data.table(pubemp.full[which(pubemp.full$confed.south == 1),])
  out2 <- out2[is.na(out2$year) == 0 , list(totadmin = mean(govemp, na.rm = T),
                                            police = mean(police, na.rm = T)), by = list(pc, year)]

# Plot Trends
# Aggregate Trends

  pubemptrends <- ggplot(out2, aes(x = as.numeric(year), y = totadmin, group = factor(pc))) + 
    geom_line(data = out2, aes(linetype = factor(pc))) + scale_x_continuous(breaks = seq(1956, 1988, by = 2)) +
    geom_vline(xintercept = 1965, color = "red") + geom_vline(xintercept = 1969, color = "red") + ggtitle("General Control and Financial Administration Employees") + xlab("") + 
    annotate("text", x = 1961, y = 60, label = "Voting Rights Act \n (1965)", color = "red") +
    annotate("text", x = 1978, y = 47, label = "Allen v. State Board of \n Elections (1969)", color = "red") +
    ylab("Average # of General Control and Financial \n Administration Employees per County") + labs(linetype = "Preclearance") + theme_bw() + 
    theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
  
  ggsave(plot = pubemptrends, filename = "ms_figure_1_right.pdf", width = 7, height = 5)


  policetrends <- ggplot(out2, aes(x = as.numeric(year), y = police, group = factor(pc))) + 
    geom_line(data = out2, aes(linetype = factor(pc))) + scale_x_continuous(breaks = seq(1956, 1988, by = 2)) +
    geom_vline(xintercept = 1965, color = "red") + geom_vline(xintercept = 1969, color = "red") + ggtitle("Police Employees") + xlab("") + 
    annotate("text", x = 1960, y = 55, label = "Voting Rights Act \n (1965)", color = "red") +
    annotate("text", x = 1974, y = 55, label = "Allen v. State Board of \n Elections (1969)", color = "red") +
    ylab("Average # of Police Employees per County") + labs(linetype = "Preclearance") + theme_bw() + 
    theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
  
  ggsave(plot = policetrends, filename = "si_figure_8.pdf", width = 7, height = 5)

##################################################################################################################################################################################################################
# Panelmatch
##################################################################################################################################################################################################################
# Get a panel match  
  longdatpm <- as.data.frame(longdat[which(is.na(longdat$totelected) == F & longdat$confed.south == 1),])
  
  longdatpm <- longdatpm %>%
    mutate(period = case_when(
      year == 1950 ~ 1,
      year == 1960 ~ 2,
      year == 1970 ~ 3,
      year == 1980 ~ 4,
      year == 1990 ~ 5,
      TRUE ~ 6))
  
  longdatpm$period <- as.integer(longdatpm$period)
  longdatpm$treated <- longdatpm$preclearance * longdatpm$vra
  longdatpm$unitid <- as.numeric(longdatpm$fips)
  
  pdf("si_figure_10.pdf")
  DisplayTreatment(unit.id = "unitid",
                   time.id = "period", legend.position = "none",
                   xlab = "Year", ylab = "County", x.angle = 0,
                   treatment = "treated", data = longdatpm, y.size = 0.5)
  dev.off()
  
##################################################################################################################################################################################################################
# Jacknife Out Individual States
##################################################################################################################################################################################################################
# Remove States Individually and Recalculate Results
  statesvec <- unique(longdat$state_name[longdat$confed.south == 1])
  
  jkstate <- data.frame(estimate = numeric(0), clusterse = numeric(0), lbound95 = numeric(0), ubound95 = numeric(0), lbound90 = numeric(0),
                        ubound90 = numeric(0), vars = character(0), label = character(0), statedropped = character(0))  
  
  # Loop over states
  for(i in 1:length(statesvec)){
    tempdat <- longdat[which(!(longdat$state_name %in% statesvec[i])),]
    
    # Base Results  
    baseres.jack <- felm(totelected ~ preclearance:vra | factor(fips) + factor(year) | 0 | state_name,  data = tempdat[tempdat$year %in% c(1950, 1960) & tempdat$confed.south == 1,])
    
    ptblackpop60.jack <- felm(totelected ~ preclearance:vra:pctblack60 | factor(fips) + factor(year) | 0 | state_name,  data = tempdat[tempdat$year %in% c(1950, 1960, 1970) & tempdat$confed.south == 1,])
    
    # Store Results    
    dropstateres <- as.data.frame(rbind(summary(baseres.jack)$coefficients["preclearance:vra",], summary(ptblackpop60.jack)$coefficients["preclearance:vra:pctblack60",]))
    
    
    dropstateres$lbound95 <- dropstateres$Estimate - qnorm(.975) * dropstateres$`Cluster s.e.`
    dropstateres$ubound95 <- dropstateres$Estimate + qnorm(.975) * dropstateres$`Cluster s.e.`
    
    dropstateres$lbound90 <- dropstateres$Estimate - qnorm(.95) * dropstateres$`Cluster s.e.`
    dropstateres$ubound90 <- dropstateres$Estimate + qnorm(.95) * dropstateres$`Cluster s.e.`
    
    dropstateres$vars <- c("Base", "Interaction")
    dropstateres$label <- dropstateres$vars
    dropstateres$statedropped <- statesvec[i]
    
    dropstateres <- dropstateres[, c(1, 2, 5, 6, 7, 8, 9, 10, 11)]
    colnames(dropstateres) <- colnames(jkstate)
    
    jkstate <- as.data.frame(rbind(jkstate, dropstateres))
    
    print(statesvec[i])
  }
  
  jkstate$label <- paste(jkstate$label, jkstate$statedropped, sep = " ")
  
  maindrop <- ggplot(data = jkstate[grepl("Base", jkstate$vars) == 1,], aes(group = statedropped)) + 
    geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
    geom_errorbarh(aes(y = statedropped, xmin = lbound95, xmax = ubound95, height = 0.1, alpha = 0.7)) + 
    geom_errorbarh(aes(y = statedropped, xmin = lbound90, xmax = ubound90, height = 0.1, alpha = 1)) + 
    geom_point(aes(x = estimate, y = statedropped), size = 3) + 
    ylab("") + xlab("Difference in Differences Estimate") +
    ggtitle("Jackknifed Estimates: Reduction in Total Elected Officials") + 
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
  
  ggsave(plot = maindrop, filename = "si_figure_12.pdf", width = 7, height = 5)
  
  interdrop <- ggplot(data = jkstate[grepl("Interaction", jkstate$vars) == 1,], aes(group = statedropped)) + 
    geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
    geom_errorbarh(aes(y = statedropped, xmin = lbound95, xmax = ubound95, height = 0.1, alpha = 0.7)) + 
    geom_errorbarh(aes(y = statedropped, xmin = lbound90, xmax = ubound90, height = 0.1, alpha = 1)) + 
    geom_point(aes(x = estimate, y = statedropped), size = 3) + 
    ylab("") + xlab("Triple Differences Estimate") +
    ggtitle("Jackknifed Estimates: Additional Reduction in \n Total Elected Officials for Counties with Larger Black Populations") + 
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5), legend.position = "none")  
  
  ggsave(plot = interdrop, filename = "si_figure_13.pdf", width = 7, height = 5)
  
##################################################################################################################################################################################################################
# Border Counties Results
##################################################################################################################################################################################################################
# Border Counties  
  # Load Data Identifying Border Counties  
  bcs <- read.csv("bordercounties.csv", header = T, stringsAsFactors = F)
  
  # Generate Border County Key to flag these in full data  
  bordercs <- c(paste(toupper(bcs$state), toupper(bcs$county), sep = "."), paste(toupper(bcs$border.state), toupper(bcs$border.county), sep = "."))
  longdat$bordercounty <- ifelse(paste(toupper(longdat$state_name), toupper(longdat$county_name), sep = ".") %in% bordercs, 1, 0)
  
  bordercounties.base <- felm(totelected ~ preclearance:vra | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$bordercounty == 1,])
  bordercounties.inter <- felm(totelected ~ preclearance:vra:pctblack60 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[longdat$year %in% c(1950, 1960) & longdat$bordercounty == 1,])
  
  stargazer(bordercounties.base, bordercounties.inter, digits = 2, title = "Difference-in-Differences Results for Border Counties", 
            out = "si_table_18.tex", label = "tab:bcs_res",
            table.placement = "h!", dep.var.labels = "County Elected Officials", 
            notes = c("Standard errors clustered by state."),
            covariate.labels = c("Preclearance x VRA", "Preclearance x VRA x Black Population in 1960"),
            single.row = T,
            no.space = T, font.size = "footnotesize")
  
##################################################################################################################################################################################################################
# Matching
##################################################################################################################################################################################################################
# Match Across Covariates
  matchstart <- fd[fd$confed.south == 1, c("preclearance", "totpop", "popdensity", "pctfem", "pct65plus", "pctblack1960", "fips", "county_name", "state_name")]
  
  matchstart <- matchstart[complete.cases(matchstart),] # this package can't match missing data so we'll just reduce to complete cases
  
  nearest.match <- matchit(formula = preclearance ~ I(totpop / 1000) + popdensity + pctfem + pct65plus + pctblack1960,
                           data = matchstart,
                           method = "nearest", 
                           distance = "logit",
                           discard = "control")
  
# Extract Data from Matching Process
  matchdat <- match.data(nearest.match)
  
# Matched Model
  matched.base <- felm(totelected ~ preclearance:vra | factor(fips) + factor(year) | 0 | state_name,  data = longdat[which(longdat$fips %in% matchdat$fips & longdat$year %in% c(1950, 1960)),])
  matched.inter <- felm(totelected ~ preclearance:vra:pctblack60 | factor(fips) + factor(year) | 0 | state_name,  data = longdat[which(longdat$fips %in% matchdat$fips & longdat$year %in% c(1950, 1960)),])
  
  stargazer(matched.base, matched.inter, digits = 2, title = "Difference-in-Differences Results with Matched Counties", 
            out = "si_table_19.tex", label = "tab:matching",
            table.placement = "h!", dep.var.labels = "County Elected Officials", 
            notes = c("Standard errors clustered by state."),
            covariate.labels = c("Preclearance x VRA", "Preclearance x VRA x Black Population in 1960"),
            single.row = T,
            no.space = T, font.size = "footnotesize")
  
##################################################################################################################################################################################################################
# Population Map
##################################################################################################################################################################################################################
  mapcts <- usmap::us_map(regions = "counties")
  mapcts <- as_tibble(merge(mapcts, longdat[which(longdat$year == 1990), c("fips", "pctblack60")], by = "fips", all.x = T))
  
  
  popmap60 <- plot_usmap(data = mapcts[, c("fips", "pctblack60")], values = "pctblack60", color = "grey90") + 
    scale_fill_continuous(name = "1960 County Black Population Share", high = "black", low = "grey90") + 
    theme(legend.position = "bottom") + 
    geom_polygon(data = usmapdata::us_map(regions = "states"),
                 aes(x, y, group = group), fill = NA, size = 1, color = "gray90")
  
  ggsave(plot = popmap60, filename = "ms_figure_7.pdf", width = 7, height = 5)

##################################################################################################################################################################################################################
# State Level Results
##################################################################################################################################################################################################################
# Collapse to State Level
  states <- data.table(longdat[which(longdat$confed.south == 1),])
  states <- states[ , list(totelected = mean(totelected, na.rm = T),
                           pctblack60 = mean(pctblack60, na.rm = T)), by = list(state_name, year)]
  
  states$vra <- ifelse(states$year > 1950, 1, 0)
  
  states$preclearance <- ifelse(states$state_name %in% c("Alabama", "Georgia", "Louisiana", "Mississippi", "South Carolina", "Virginia", "North Carolina") & states$year > 1950, 1, 0)
  states$preclearance2 <- ifelse(states$state_name %in% c("Alabama", "Georgia", "Louisiana", "Mississippi", "South Carolina", "Virginia") & states$year > 1950, 1, 0)
  
  statereg <- felm(totelected ~ preclearance:vra | factor(state_name) + factor(year),  data = states[states$year %in% c(1950, 1960),])
  statereg2 <- felm(totelected ~ preclearance2:vra | factor(state_name) + factor(year),  data = states[states$year %in% c(1950, 1960),])
  
  statereg.inter <- felm(totelected ~ preclearance:vra:pctblack60 | factor(state_name) + factor(year),  data = states[states$year %in% c(1950, 1960),])
  statereg.inter2 <- felm(totelected ~ preclearance2:vra:pctblack60 | factor(state_name) + factor(year),  data = states[states$year %in% c(1950, 1960),])
  
  stargazer(statereg, statereg.inter, statereg2, statereg.inter2, digits = 2, title = "State Level Difference-in-Differences Results", 
            out = "si_table_9.tex", label = "tab:statedid",
            table.placement = "H", dep.var.labels = "County Elected Officials", covariate.labels = c("Covered x VRA, NC = 1", "Covered x VRA x 1960 Black Pop., NC = 1", "Covered x VRA, NC = 0", "Covered x VRA x 1960 Black Population, NC = 0"),
            no.space = T, font.size = "footnotesize")

